#######		REPLICATION FILES
#######		Conflict versus Disaster-induced Displacement: Similar or Distinct Implications for Security?
#######		Civil Wars 
#######		Heidrun Bohnet, Fabien Cottier and Simon Hug
#######		Replication: Results of the empirical analyses shown in the appendix
#######		This version: September 2021


#### setup

# set working directory and load required packages

rm(list = ls())
setwd("path/to/Working/directory")
library("foreign")
library("lmtest")
library("sandwich")
library("survival")
library("stargazer")
library("plotrix")
library("arm")
library("plyr")

# function to extract quantities of interest
source("Rscripts/qi_linear.R")
source("Rscripts/qi_interaction.R")
source("Rscripts/qi_interaction_noidp.R")

# load data
cvdm_data <- read.csv("cvdm_data.csv")
names(cvdm_data)

# create dataset for full (Model 1: 1991-2011) and restricted sample (Models M2-M4) (only period 2008-2011)
cvdm_full_df <- cvdm_data
cvdm_restricted_df <- cvdm_data[cvdm_data$year>2007,]

# creation of interaction terms (Model M4)
# full overlap admin unit
cvdm_restricted_df$iterm1_pw <- cvdm_restricted_df$DIDP_1000pw*cvdm_restricted_df$idp_ra
cvdm_restricted_df$iterm2_pw <- cvdm_restricted_df$floodpw*cvdm_restricted_df$idp_ra
cvdm_restricted_df$iterm1_pw_1tl <- cvdm_restricted_df$DIDP_1000pw_1tl*cvdm_restricted_df$idp_ra
cvdm_restricted_df$iterm2_pw_1tl <- cvdm_restricted_df$floodpw_1tl*cvdm_restricted_df$idp_ra
# partial overlap admin unit
cvdm_restricted_df$iterm1_npw <- cvdm_restricted_df$DIDP_1000npw*cvdm_restricted_df$idp_ra
cvdm_restricted_df$iterm2_npw <- cvdm_restricted_df$floodnpw*cvdm_restricted_df$idp_ra
cvdm_restricted_df$iterm1_npw_1tl <- cvdm_restricted_df$DIDP_1000npw_1tl*cvdm_restricted_df$idp_ra
cvdm_restricted_df$iterm2_npw_1tl <- cvdm_restricted_df$floodnpw_1tl*cvdm_restricted_df$idp_ra
# near admin unit (with 25 km)
cvdm_restricted_df$iterm1_near <- cvdm_restricted_df$DIDP_1000near*cvdm_restricted_df$idp_ra
cvdm_restricted_df$iterm2_near <- cvdm_restricted_df$floodnear*cvdm_restricted_df$idp_ra
cvdm_restricted_df$iterm1_near_1tl <- cvdm_restricted_df$DIDP_1000near_1tl*cvdm_restricted_df$idp_ra
cvdm_restricted_df$iterm2_near_1tl <- cvdm_restricted_df$floodnear_1tl*cvdm_restricted_df$idp_ra



#### Section 2: summary statistics and additional analysis (Models 1 and 4 in the main text)

# obtain estimates for Models 1 and 4 in the main text
model1 <- glm(sc_inc_all ~ DIDP_1000pw + floodpw + # full overlap (t0)
    DIDP_1000npw + floodnpw + # part. overlap (t0)
    DIDP_1000near + floodnear + # near (t0)
    DIDP_1000pw_1tl + floodpw_1tl + # full overlap (t-1)
    DIDP_1000npw_1tl + floodnpw_1tl + # part. overlap (t-1)
    DIDP_1000near_1tl + floodnear_1tl + # near (t-1)
    flood_pexp + flood_pexp_sq + # count past events (all events)
    agdppc_ln + iapop_ln + sc_inc_all_1tl  + sc_cl_inc_1tl + # CV admin 
    xpolity_1tl + anocracy_1tl + crgdppc_ln_1tl + cgdpg_1tl + cpop_ln_1tl + cpopg_1tl + incidencev412 + # CV country
    peaceyears_sc + peaceyears_sc_sq + peaceyears_sc_cb, # peaceyears variables,
  family = binomial(link = "logit"),
  data = cvdm_full_df)

model4 <- bayesglm(sc_inc_all ~ DIDP_1000pw + iterm1_pw + floodpw + iterm2_pw + #full overlap (t0) 
    DIDP_1000npw + iterm1_npw + floodnpw + iterm2_npw + #part. overlap (t0)
    DIDP_1000near + iterm1_near + floodnear + iterm2_near + #near (t0)
    DIDP_1000pw_1tl + iterm1_pw_1tl + floodpw_1tl + iterm2_pw_1tl +#full overlap (t-1)
    DIDP_1000npw_1tl + iterm1_npw_1tl + floodnpw_1tl + iterm2_npw_1tl + #part. overlap (t-1)
    DIDP_1000near_1tl + iterm1_near_1tl + floodnear_1tl + iterm2_near_1tl + #near (t-1)
    idp_ra + # conflict IDPs
    flood_pexp + flood_pexp_sq + # count past events (all events)
    agdppc_ln + iapop_ln + sc_inc_all_1tl  + sc_cl_inc_1tl + # CV admin 
    xpolity_1tl + anocracy_1tl + crgdppc_ln_1tl + cgdpg_1tl + cpop_ln_1tl + cpopg_1tl + incidencev412 + # CV country
    peaceyears_sc + peaceyears_sc_sq + peaceyears_sc_cb, # peaceyears variables,
  family = binomial(link = "logit"),
  data = cvdm_restricted_df)

## Table 1: summary statistics

model1_df<-expand.model.frame(model1,c("pe_inc_all","vi_inc_all","sc_count_all",
  "pe_inc_all_1tl","vi_inc_all_1tl","pe_cl_inc_1tl","vi_cl_inc_1tl",
  "ns_inc_1tl","ns_cl_inc_1tl","sb_inc_1tl","sb_cl_inc_1tl","ns_inc","sb_inc"),na.expand = T)
stargazer(model1_df[,c(1,32,30:31,41:42,2,4,6,3,5,7,14,16:17,19,35:36,38,40,20:27)],
  covariate.labels=c("social conflict incidence","social conflict count",
  "Peaceful conflict incidence","Violence conflict incidence","Non-state conflict incidence","Civil war incidence",
  "Disaster displ. (full overlap)","Disaster displ. (partial overlap)","Disaster displ. (near)",
  "Flood (full overlap)","Flood (partial overlap)","Flood (near)",
  "Past flood count",
  "Admin gdppc (ln)","Admin pop (ln)",
  "Spatial lag (social conflict)","Spatial lag (peaceful conflict)","Spatial lag (violent conflict)",
  "Spatial lag (non-state conflict)","Spatial lag (civil war)",
  "Xpolity (lag)","Anocracy (lag)","Country gdppc (ln lag)","Country gdp growth (lag)",
  "Country pop (ln lag)","Country pop growth (lag)","civil war","Peace years (social conflict)"),
  summary.stat = c("N","mean","median","sd","min","max"), type = "text")


## Figure A.1: : Average simulated probabilities of social conflict: past exposure to floods

# obtain quantities of interest
model1.flood_pexp<-model.matrix(model1)
coef_m1 <- coef(model1)
vcov_m1 <- vcovHC(model1)

# simulate for each number of past floods (min = 0, max = 19)
set.seed(1234)
simb_m1 <- mvrnorm(1000,coef_m1,vcov_m1)
seqsim <- seq(from=0,to=19,length.out=11)
yhat_m1 <- matrix(NA,1000,length(seqsim))
for (i in 1:length(seqsim)){
  model1.flood_pexp[,"flood_pexp"] <- seqsim[i]
  model1.flood_pexp[,"flood_pexp_sq"] <- (seqsim[i])^2
  for (j in 1:1000){
    yhat_m1[j,i] <- mean(pnorm(simb_m1[j,]%*%t(model1.flood_pexp)))
  }
} 
yhat_m1_quantile=apply(yhat_m1,2,quantile,c(0.025,0.5,0.975))

# create plot
plot(seqsim,yhat_m1_quantile[3,],ylim=range(0,0.2),type="n",main="",xlab="Number of past flood events",ylab="P (social conflict)")
for (j in 1:1000){
  points(seqsim,yhat_m1[j,], col=rgb(0,130,0,5,maxColorValue=255), type="l")
}
abline(h=0,lty=5)
lines(seqsim,yhat_m1_quantile[2,],lty=1,col=rgb(255,0,0,100,maxColorValue=255))
lines(seqsim,yhat_m1_quantile[3,],lty=2,col=rgb(0,0,130,100,maxColorValue=255))
lines(seqsim,yhat_m1_quantile[1,],lty=2,col=rgb(0,0,130,100,maxColorValue=255))
rug(jitter(model.matrix(model1)[,"flood_pexp"], amount = 0.5), side = 3,col="grey")


## Figure A.2: : Average predictive differences in conflict probabilities in the absence of conflict IDPs (Model 4)

# obtain quantities of interest
apd_matrix_m4_noidp_0tl <- qi_inter_noidp(model4,lag=F)
apd_matrix_m4_noidp_1tl <- qi_inter_noidp(model4,lag=T)

# create plot
yat=rev(c(seq(0,2/3,by = 1/3),seq(0,2/3,by = 1/3)+(2/3+.5),seq(0,2/3,by = 1/3)+(4/3+1)))
limx=c(-.3,.3)
limx_negbin=c(-1,1)
par(ps=15,las=1,mar=c(0,0,0,2),oma=c(4, 12, 3, 0),mfcol=c(1,2))
# plot 1
plot.new();plot.window(xlim = limx,ylim=range(yat));abline(v=0,col="grey")
plotCI(apd_matrix_m4_noidp_1tl[,1],yat,
        li=apd_matrix_m4_noidp_1tl[,3],ui=apd_matrix_m4_noidp_1tl[,6],err='x',
        gap=0,xlab="",bty="n",
        ylab="",yaxt="n",pch=16,add=T);axis(1,cex.axis=.6)
axis(side=2, at=yat, labels=c(
    "Fully affected, Floods","Fully affected, Displ.",expression(italic("Difference")),
    "Partially affected, Floods","Partially affected, Displ.",expression(italic("Difference")),
    "Proximate, Floods","Proximate, Displ.",expression(italic("Difference"))
    ),cex.axis=.6,tick=FALSE)
mtext(expression(paste(Delta," P (Social conflict)")), 1, 2.5, cex=.6, outer=FALSE)
mtext("Lagged response (T-1)", 3, 1, cex=.6,at=0, outer=FALSE)
# plot 2
plot.new();plot.window(xlim = limx,ylim=range(yat));abline(v=0,col="grey")
plotCI(apd_matrix_m4_noidp_0tl[,1],yat,
        li=apd_matrix_m4_noidp_0tl[,3],ui=apd_matrix_m4_noidp_0tl[,6],err='x',
        gap=0,xlab="",bty="n",
        ylab="",yaxt="n",pch=16,add=T);axis(1,cex.axis=.6)
mtext(expression(paste(Delta," P (Social conflict)")), 1, 2.5, cex=.6, outer=FALSE)
mtext("Immediate response (T0)", 3, 1, cex=.6,at=0, outer=FALSE)





#### sensitivity analyses: Model 1

# adjust graph plot window limits (Figures A.8--11)
yat=rev(c(seq(0,2/3,by = 1/3),seq(0,2/3,by = 1/3)+(2/3+.5),seq(0,2/3,by = 1/3)+(4/3+1)))
limx=c(-.3,.3)
limx_negbin=c(-1,1)

## Model s1.m1: all peaceful events
s1.m1 <- glm(pe_inc_all ~ DIDP_1000pw + floodpw + # full overlap (t0)
    DIDP_1000npw + floodnpw + # part. overlap (t0)
    DIDP_1000near + floodnear + # near (t0)
    DIDP_1000pw_1tl + floodpw_1tl + # full overlap (t-1)
    DIDP_1000npw_1tl + floodnpw_1tl + # part. overlap (t-1)
    DIDP_1000near_1tl + floodnear_1tl + # near (t-1)
    flood_pexp + flood_pexp_sq + # count past events (all events)
    agdppc_ln + iapop_ln + pe_inc_all_1tl  + pe_cl_inc_1tl + # CV admin 
    xpolity_1tl + anocracy_1tl + crgdppc_ln_1tl + cgdpg_1tl + cpop_ln_1tl + cpopg_1tl + incidencev412 + # CV country
    peaceyears_pe + peaceyears_pe_sq + peaceyears_pe_cb, # time dep
  family = binomial(link="logit"),
  data = cvdm_full_df)
coeftest(s1.m1,vcovHC(s1.m1,"HC0"))
# N obs
nobs(s1.m1)
# Loglik
logLik(s1.m1)
# aic
AIC(s1.m1)

# Figure A.3:  average predictive differences
apd_matrix_m1S1_0tl<-qi_linear(s1.m1,lag=F)
apd_matrix_m1S1_1tl<-qi_linear(s1.m1,lag=T)
par(ps=15,las=1,mar=c(0,0,0,2),oma=c(4, 12, 3, 0),mfcol=c(1,2))
# plot 1
plot.new();plot.window(xlim = limx,ylim=range(yat));abline(v=0,col="grey")
plotCI(apd_matrix_m1S1_1tl[,1],yat,
    li=apd_matrix_m1S1_1tl[,3],ui=apd_matrix_m1S1_1tl[,6],err='x',
    gap=0,xlab="",bty="n",
    ylab="",yaxt="n",pch=16,add=T);axis(1,cex.axis=.6)
axis(side=2, at=yat, labels=c(
    "Fully affected, Floods","Fully affected, Displ.",expression(italic("Difference")),
    "Partially affected, Floods","Partially affected, Displ.",expression(italic("Difference")),
    "Proximate, Floods","Proximate, Displ.",expression(italic("Difference"))
    ),cex.axis=.6,tick=FALSE)
mtext(expression(paste(Delta," P (Peaceful social conflict)")), 1, 2.5, cex=.6, outer=FALSE)
mtext("Lagged response (T-1)", 3, 1, cex=.6,at=0, outer=FALSE)
# plot 2
plot.new();plot.window(xlim = limx,ylim=range(yat));abline(v=0,col="grey")
plotCI(apd_matrix_m1S1_0tl[,1],yat,
    li=apd_matrix_m1S1_0tl[,3],ui=apd_matrix_m1S1_0tl[,6],err='x',
    gap=0,xlab="",bty="n",
    ylab="",yaxt="n",pch=16,add=T);axis(1,cex.axis=.6)
mtext(expression(paste(Delta," P (Peaceful social conflict)")), 1, 2.5, cex=.6, outer=FALSE)
mtext("Immediate response (T0)", 3, 1, cex=.6,at=0, outer=FALSE)


## Model s2.m1: all violent events
s2.m1 <- glm(vi_inc_all ~ DIDP_1000pw + floodpw + # full overlap (t0)
    DIDP_1000npw + floodnpw + # part. overlap (t0)
    DIDP_1000near + floodnear + # near (t0)
    DIDP_1000pw_1tl + floodpw_1tl + # full overlap (t-1)
    DIDP_1000npw_1tl + floodnpw_1tl + # part. overlap (t-1)
    DIDP_1000near_1tl + floodnear_1tl + # near (t-1)
    flood_pexp + flood_pexp_sq + # count past events (all events)
    agdppc_ln + iapop_ln + vi_inc_all_1tl  + vi_cl_inc_1tl + # CV admin 
    xpolity_1tl + anocracy_1tl + crgdppc_ln_1tl + cgdpg_1tl + cpop_ln_1tl + cpopg_1tl + incidencev412 + # CV country
    peaceyears_vi + peaceyears_vi_sq + peaceyears_vi_cb, # time dep
  family = binomial(link="logit"),
  data = cvdm_full_df)
coeftest(s2.m1,vcovHC(s2.m1,"HC0"))
# N obs
nobs(s2.m1)
# Loglik
logLik(s2.m1)
# aic
AIC(s2.m1)

# Figure A.4:  average predictive differences
apd_matrix_m1S2_0tl<-qi_linear(s2.m1,lag=F)
apd_matrix_m1S2_1tl<-qi_linear(s2.m1,lag=T)
par(ps=15,las=1,mar=c(0,0,0,2),oma=c(4, 12, 3, 0),mfcol=c(1,2))
# plot 1
plot.new();plot.window(xlim = limx,ylim=range(yat));abline(v=0,col="grey")
plotCI(apd_matrix_m1S2_1tl[,1],yat,
    li=apd_matrix_m1S2_1tl[,3],ui=apd_matrix_m1S2_1tl[,6],err='x',
    gap=0,xlab="",bty="n",
    ylab="",yaxt="n",pch=16,add=T);axis(1,cex.axis=.6)
axis(side=2, at=yat, labels=c(
    "Fully affected, Floods","Fully affected, Displ.",expression(italic("Difference")),
    "Partially affected, Floods","Partially affected, Displ.",expression(italic("Difference")),
    "Proximate, Floods","Proximate, Displ.",expression(italic("Difference"))
    ),cex.axis=.6,tick=FALSE)
mtext(expression(paste(Delta," P (Violent social conflict)")), 1, 2.5, cex=.6, outer=FALSE)
mtext("Lagged response (T-1)", 3, 1, cex=.6,at=0, outer=FALSE)
# plot 2
plot.new();plot.window(xlim = limx,ylim=range(yat));abline(v=0,col="grey")
plotCI(apd_matrix_m1S2_0tl[,1],yat,
    li=apd_matrix_m1S2_0tl[,3],ui=apd_matrix_m1S2_0tl[,6],err='x',
    gap=0,xlab="",bty="n",
    ylab="",yaxt="n",pch=16,add=T);axis(1,cex.axis=.6)
mtext(expression(paste(Delta," P (Violent social conflict)")), 1, 2.5, cex=.6, outer=FALSE)
mtext("Immediate response (T0)", 3, 1, cex=.6,at=0, outer=FALSE)


## Model s3.m1: negative binomial
s3.m1 <- glm.nb(sc_count_all ~ DIDP_1000pw + floodpw + # full overlap (t0)
    DIDP_1000npw + floodnpw + # part. overlap (t0)
    DIDP_1000near + floodnear + # near (t0)
    DIDP_1000pw_1tl + floodpw_1tl + # full overlap (t-1)
    DIDP_1000npw_1tl + floodnpw_1tl + # part. overlap (t-1)
    DIDP_1000near_1tl + floodnear_1tl + # near (t-1)
    flood_pexp + flood_pexp_sq + # count past events (all events)
   	agdppc_ln + iapop_ln + sc_count_all_1tl  + sc_cl_count_1tl + # CV admin 
    xpolity_1tl + anocracy_1tl + crgdppc_ln_1tl + cgdpg_1tl + cpop_ln_1tl + cpopg_1tl + incidencev412 + # CV country
    peaceyears_sc + peaceyears_sc_sq + peaceyears_sc_cb, # time dep
  data = cvdm_full_df)
coeftest(s3.m1,vcovHC(s3.m1,"HC0"))
# N obs
nobs(s3.m1)
# Dispersion parameter (theta) + associated SE
s3.m1$theta; s3.m1$SE.theta
# Loglik
logLik(s3.m1)
# aic
AIC(s3.m1)

# Figure A.5:  average predictive differences
apd_matrix_m1S3_0tl<-qi_linear(s3.m1,lag=F)
apd_matrix_m1S3_1tl<-qi_linear(s3.m1,lag=T)
par(ps=15,las=1,mar=c(0,0,0,2),oma=c(4, 12, 3, 0),mfcol=c(1,2))
# plot 1
plot.new();plot.window(xlim = limx_negbin,ylim=range(yat));abline(v=0,col="grey")
plotCI(apd_matrix_m1S3_1tl[,1],yat,
    li=apd_matrix_m1S3_1tl[,3],ui=apd_matrix_m1S3_1tl[,6],err='x',
    gap=0,xlab="",bty="n",
    ylab="",yaxt="n",pch=16,add=T);axis(1,cex.axis=.6)
axis(side=2, at=yat, labels=c(
    "Fully affected, Floods","Fully affected, Displ.",expression(italic("Difference")),
    "Partially affected, Floods","Partially affected, Displ.",expression(italic("Difference")),
    "Proximate, Floods","Proximate, Displ.",expression(italic("Difference"))
    ),cex.axis=.6,tick=FALSE)
mtext(expression(paste(Delta," N (Social conflict)")), 1, 2.5, cex=.6, outer=FALSE)
mtext("Lagged response (T-1)", 3, 1, cex=.6,at=0, outer=FALSE)
# plot 2
plot.new();plot.window(xlim = limx_negbin,ylim=range(yat));abline(v=0,col="grey")
plotCI(apd_matrix_m1S3_0tl[,1],yat,
    li=apd_matrix_m1S3_0tl[,3],ui=apd_matrix_m1S3_0tl[,6],err='x',
    gap=0,xlab="",bty="n",
    ylab="",yaxt="n",pch=16,add=T);axis(1,cex.axis=.6)
mtext(expression(paste(Delta," N (Social conflict)")), 1, 2.5, cex=.6, outer=FALSE)
mtext("Immediate response (T0)", 3, 1, cex=.6,at=0, outer=FALSE)


## Model s4.m1: conditional fixed-effect logit
s4.m1 <- clogit(sc_inc_all ~ DIDP_1000pw + floodpw + # full overlap (t0) 
    DIDP_1000npw + floodnpw + # part. overlap (t0) 
    DIDP_1000near + floodnear + # near (t0) 
    DIDP_1000pw_1tl + floodpw_1tl + # full overlap (t-1) 
    DIDP_1000npw_1tl + floodnpw_1tl + # part. overlap (t-1) 
    DIDP_1000near_1tl + floodnear_1tl + # near (t-1) 
    flood_pexp + flood_pexp_sq + # count past events (all events)
    agdppc_ln + iapop_ln + sc_inc_all_1tl  + sc_cl_inc_1tl + # CV admin 
    xpolity_1tl + anocracy_1tl + crgdppc_ln_1tl + cgdpg_1tl + cpop_ln_1tl + cpopg_1tl + incidencev412 + # CV country
    peaceyears_sc + peaceyears_sc_sq + peaceyears_sc_cb + # time dep
   strata(admin1_code),
   data = cvdm_full_df)
summary(s4.m1)
# N obs
nrow(model.matrix(s4.m1))
# Loglik
logLik(s4.m1)
# aic
AIC(s4.m1)


## Model s5.m1 7: UCDP non-state conflict
s5.m1 <- glm(ns_inc ~ DIDP_1000pw + floodpw + # full overlap (t0) 
    DIDP_1000npw + floodnpw + # part. overlap (t0) 
    DIDP_1000near + floodnear + # near (t0) 
    DIDP_1000pw_1tl + floodpw_1tl + # full overlap (t-1) 
    DIDP_1000npw_1tl + floodnpw_1tl + # part. overlap (t-1) 
    DIDP_1000near_1tl + floodnear_1tl + # near (t-1) 
    flood_pexp + flood_pexp_sq + # count past events (all events)
    agdppc_ln + iapop_ln + ns_inc_1tl  + ns_cl_inc_1tl + # CV admin 
    xpolity_1tl + anocracy_1tl + crgdppc_ln_1tl + cgdpg_1tl + cpop_ln_1tl + cpopg_1tl + incidencev412 + # CV country
    peaceyears_ns + peaceyears_ns_sq + peaceyears_ns_cb, # time dep
  family = binomial(link="logit"),
  data = cvdm_full_df)
coeftest(s5.m1,vcovHC(s5.m1,"HC0"))
# N obs
nobs(s5.m1)
# Loglik
logLik(s5.m1)
# aic
AIC(s5.m1)

# Figure A.6: average predictive differences
apd_matrix_m1S5_0tl <- qi_linear(s5.m1,lag=F)
apd_matrix_m1S5_1tl <- qi_linear(s5.m1,lag=T)
par(ps=15,las=1,mar=c(0,0,0,2),oma=c(4, 12, 3, 0),mfcol=c(1,2))
# plot 1
plot.new();plot.window(xlim = limx,ylim=range(yat));abline(v=0,col="grey")
plotCI(apd_matrix_m1S5_1tl[,1],yat,
    li=apd_matrix_m1S5_1tl[,3],ui=apd_matrix_m1S5_1tl[,6],err='x',
    gap=0,xlab="",bty="n",
    ylab="",yaxt="n",pch=16,add=T);axis(1,cex.axis=.6)
axis(side=2, at=yat, labels=c(
    "Fully affected, Floods","Fully affected, Displ.",expression(italic("Difference")),
    "Partially affected, Floods","Partially affected, Displ.",expression(italic("Difference")),
    "Proximate, Floods","Proximate, Displ.",expression(italic("Difference"))
    ),cex.axis=.6,tick=FALSE)
mtext(expression(paste(Delta," P (Social conflict)")), 1, 2.5, cex=.6, outer=FALSE)
mtext("Lagged response (T-1)", 3, 1, cex=.6,at=0, outer=FALSE)
# plot 2
plot.new();plot.window(xlim = limx,ylim=range(yat));abline(v=0,col="grey")
plotCI(apd_matrix_m1S5_0tl[,1],yat,
    li=apd_matrix_m1S5_0tl[,3],ui=apd_matrix_m1S5_0tl[,6],err='x',
    gap=0,xlab="",bty="n",
    ylab="",yaxt="n",pch=16,add=T);axis(1,cex.axis=.6)
mtext(expression(paste(Delta," P (Social conflict)")), 1, 2.5, cex=.6, outer=FALSE)
mtext("Immediate response (T0)", 3, 1, cex=.6,at=0, outer=FALSE)


## Model s6.m1: UCDP civil war conflict
s6.m1 <- glm(sb_inc ~ DIDP_1000pw + floodpw + # full overlap (t0) 
    DIDP_1000npw + floodnpw + # part. overlap (t0) 
    DIDP_1000near + floodnear + # near (t0) 
    DIDP_1000pw_1tl + floodpw_1tl + # full overlap (t-1) 
    DIDP_1000npw_1tl + floodnpw_1tl + # part. overlap (t-1) 
    DIDP_1000near_1tl + floodnear_1tl + # near (t-1) 
    flood_pexp + flood_pexp_sq + # count past events (all events)
    agdppc_ln + iapop_ln + sb_inc_1tl  + sb_cl_inc_1tl + # CV admin 
    xpolity_1tl + anocracy_1tl + crgdppc_ln_1tl + cgdpg_1tl + cpop_ln_1tl + cpopg_1tl + incidencev412 + # CV country
    peaceyears_sb + peaceyears_sb_sq + peaceyears_sb_cb, # time dep
  family = binomial(link="logit"),
  data = cvdm_full_df)
coeftest(s6.m1,vcovHC(s6.m1,"HC0"))
# N obs
nobs(s6.m1)
# Loglik
logLik(s6.m1)
# aic
AIC(s6.m1)

# Figure A.7: average predictive differences
apd_matrix_m1S6_0tl<-qi_linear(s6.m1,lag=F)
apd_matrix_m1S6_1tl<-qi_linear(s6.m1,lag=T)
par(ps=15,las=1,mar=c(0,0,0,2),oma=c(4, 12, 3, 0),mfcol=c(1,2))
# plot 1
plot.new();plot.window(xlim = limx,ylim=range(yat));abline(v=0,col="grey")
plotCI(apd_matrix_m1S6_1tl[,1],yat,
    li=apd_matrix_m1S6_1tl[,3],ui=apd_matrix_m1S6_1tl[,6],err='x',
    gap=0,xlab="",bty="n",
    ylab="",yaxt="n",pch=16,add=T);axis(1,cex.axis=.6)
axis(side=2, at=yat, labels=c(
    "Fully affected, Floods","Fully affected, Displ.",expression(italic("Difference")),
    "Partially affected, Floods","Partially affected, Displ.",expression(italic("Difference")),
    "Proximate, Floods","Proximate, Displ.",expression(italic("Difference"))
    ),cex.axis=.6,tick=FALSE)
mtext(expression(paste(Delta," P (Civil war)")), 1, 2.5, cex=.6, outer=FALSE)
mtext("Lagged response (T-1)", 3, 1, cex=.6,at=0, outer=FALSE)
# plot 2
plot.new();plot.window(xlim = limx,ylim=range(yat));abline(v=0,col="grey")
plotCI(apd_matrix_m1S6_0tl[,1],yat,
    li=apd_matrix_m1S6_0tl[,3],ui=apd_matrix_m1S6_0tl[,6],err='x',
    gap=0,xlab="",bty="n",
    ylab="",yaxt="n",pch=16,add=T);axis(1,cex.axis=.6)
mtext(expression(paste(Delta," P (Civil war)")), 1, 2.5, cex=.6, outer=FALSE)
mtext("Immediate response (T0)", 3, 1, cex=.6,at=0, outer=FALSE)



#### sensitivity analyses: Model 4

# adjust graph plot window limits (Figures A.8--11)
limx <- c(-1,1)
limx_negbin <- c(-5,5)


## Model s1.m4: peaceful events
s1.m4 <- bayesglm(pe_inc_all ~ DIDP_1000pw + iterm1_pw + floodpw + iterm2_pw +# full overlap (t0) 
    DIDP_1000npw + iterm1_npw + floodnpw + iterm2_npw + # part. overlap (t0) 
    DIDP_1000near + iterm1_near + floodnear + iterm2_near + # near (t0) 
    DIDP_1000pw_1tl + iterm1_pw_1tl + floodpw_1tl + iterm2_pw_1tl +# full overlap (t-1)
    DIDP_1000npw_1tl + iterm1_npw_1tl + floodnpw_1tl + iterm2_npw_1tl + # part. overlap (t-1)
    DIDP_1000near_1tl + iterm1_near_1tl + floodnear_1tl + iterm2_near_1tl + # near (t-1)
    idp_ra + # conflict displacement variable
    flood_pexp + flood_pexp_sq + # count past events (all events)
    agdppc_ln + iapop_ln + pe_inc_all_1tl  + pe_cl_inc_1tl + # CV admin 
    xpolity_1tl + anocracy_1tl + crgdppc_ln_1tl + cgdpg_1tl + cpop_ln_1tl + cpopg_1tl + incidencev412 + # CV country
    peaceyears_pe + peaceyears_pe_sq + peaceyears_pe_cb, # time dep
  family = binomial(link="logit"),
  data = cvdm_restricted_df)
coeftest(s1.m4,vcovHC(s1.m4,"HC0"))
# N obs
nobs(s1.m4)
# Loglik
logLik(s1.m4)
# aic
AIC(s1.m4)

# Figure A.8: average predictive differences
apd_matrix_m4S1_0tl<-qi_inter(s1.m4,lag=F)
apd_matrix_m4S1_1tl<-qi_inter(s1.m4,lag=T)
par(ps=15,las=1,mar=c(0,0,0,2),oma=c(4, 12, 3, 0),mfcol=c(1,2))
# plot 1
plot.new();plot.window(xlim = limx,ylim=range(yat));abline(v=0,col="grey")
plotCI(apd_matrix_m4S1_1tl[,1],yat,
    li=apd_matrix_m4S1_1tl[,3],ui=apd_matrix_m4S1_1tl[,6],err='x',
    gap=0,xlab="",bty="n",
    ylab="",yaxt="n",pch=16,add=T);axis(1,cex.axis=.6)
axis(side=2, at=yat, labels=c(
    "Fully affected, Floods","Fully affected, Displ.",expression(italic("Difference")),
    "Partially affected, Floods","Partially affected, Displ.",expression(italic("Difference")),
    "Proximate, Floods","Proximate, Displ.",expression(italic("Difference"))
    ),cex.axis=.6,tick=FALSE)
mtext(expression(paste(Delta," P (Peaceful social conflict)")), 1, 2.5, cex=.6, outer=FALSE)
mtext("Lagged response (T-1)", 3, 1, cex=.6,at=0, outer=FALSE)
# plot 2
plot.new();plot.window(xlim = limx,ylim=range(yat));abline(v=0,col="grey")
plotCI(apd_matrix_m4S1_0tl[,1],yat,
    li=apd_matrix_m4S1_0tl[,3],ui=apd_matrix_m4S1_0tl[,6],err='x',
    gap=0,xlab="",bty="n",
    ylab="",yaxt="n",pch=16,add=T);axis(1,cex.axis=.6)
mtext(expression(paste(Delta," P (Peaceful social conflict)")), 1, 2.5, cex=.6, outer=FALSE)
mtext("Immediate response (T0)", 3, 1, cex=.6,at=0, outer=FALSE)


## Model s2.m4: violent events
s2.m4 <- bayesglm(vi_inc_all ~ DIDP_1000pw + iterm1_pw + floodpw + iterm2_pw +# full overlap (t0) 
    DIDP_1000npw + iterm1_npw + floodnpw + iterm2_npw + # part. overlap (t0) 
    DIDP_1000near + iterm1_near + floodnear + iterm2_near + # near (t0) 
    DIDP_1000pw_1tl + iterm1_pw_1tl + floodpw_1tl + iterm2_pw_1tl +# full overlap (t-1)
    DIDP_1000npw_1tl + iterm1_npw_1tl + floodnpw_1tl + iterm2_npw_1tl + # part. overlap (t-1)
    DIDP_1000near_1tl + iterm1_near_1tl + floodnear_1tl + iterm2_near_1tl + # near (t-1)
    idp_ra + # conflict displacement variable
    flood_pexp + flood_pexp_sq + # count past events (all events)
    agdppc_ln + iapop_ln + vi_inc_all_1tl  + vi_cl_inc_1tl + # CV admin 
    xpolity_1tl + anocracy_1tl + crgdppc_ln_1tl + cgdpg_1tl + cpop_ln_1tl + cpopg_1tl + incidencev412 + # CV country
    peaceyears_vi + peaceyears_vi_sq + peaceyears_vi_cb, # time dep
  family = binomial(link="logit"),
  data = cvdm_restricted_df)
coeftest(s2.m4,vcovHC(s2.m4,"HC0"))
# N obs
nobs(s2.m4)
# Loglik
logLik(s2.m4)
# aic
AIC(s2.m4)

# Figure A.9: average predictive differences
apd_matrix_m4S2_0tl<-qi_inter(s2.m4,lag=F)
apd_matrix_m4S2_1tl<-qi_inter(s2.m4,lag=T)
par(ps=15,las=1,mar=c(0,0,0,2),oma=c(4, 12, 3, 0),mfcol=c(1,2))
# plot 1
plot.new();plot.window(xlim = limx,ylim=range(yat));abline(v=0,col="grey")
plotCI(apd_matrix_m4S2_1tl[,1],yat,
    li=apd_matrix_m4S2_1tl[,3],ui=apd_matrix_m4S2_1tl[,6],err='x',
    gap=0,xlab="",bty="n",
    ylab="",yaxt="n",pch=16,add=T);axis(1,cex.axis=.6)
axis(side=2, at=yat, labels=c(
    "Fully affected, Floods","Fully affected, Displ.",expression(italic("Difference")),
    "Partially affected, Floods","Partially affected, Displ.",expression(italic("Difference")),
    "Proximate, Floods","Proximate, Displ.",expression(italic("Difference"))
    ),cex.axis=.6,tick=FALSE)
mtext(expression(paste(Delta," P (Violent social conflict)")), 1, 2.5, cex=.6, outer=FALSE)
mtext("Lagged response (T-1)", 3, 1, cex=.6,at=0, outer=FALSE)
# plot 2
plot.new();plot.window(xlim = limx,ylim=range(yat));abline(v=0,col="grey")
plotCI(apd_matrix_m4S2_0tl[,1],yat,
    li=apd_matrix_m4S2_0tl[,3],ui=apd_matrix_m4S2_0tl[,6],err='x',
    gap=0,xlab="",bty="n",
    ylab="",yaxt="n",pch=16,add=T);axis(1,cex.axis=.6)
mtext(expression(paste(Delta," P (Violent social conflict)")), 1, 2.5, cex=.6, outer=FALSE)
mtext("Immediate response (T0)", 3, 1, cex=.6,at=0, outer=FALSE)


## Model s3.m4: negative binomial
s3.m4 <- glm.nb(sc_count_all ~ DIDP_1000pw + iterm1_pw + floodpw + iterm2_pw +# full overlap (t0) 
    DIDP_1000npw + iterm1_npw + floodnpw + iterm2_npw + # part. overlap (t0) 
    DIDP_1000near + iterm1_near + floodnear + iterm2_near + # near (t0) 
    DIDP_1000pw_1tl + iterm1_pw_1tl + floodpw_1tl + iterm2_pw_1tl +# full overlap (t-1)
    DIDP_1000npw_1tl + iterm1_npw_1tl + floodnpw_1tl + iterm2_npw_1tl + # part. overlap (t-1)
    DIDP_1000near_1tl + iterm1_near_1tl + floodnear_1tl + iterm2_near_1tl + # near (t-1)
    idp_ra + # conflict displacement variable
    flood_pexp + flood_pexp_sq + # count past events (all events)
   	agdppc_ln + iapop_ln + sc_count_all_1tl  + sc_cl_count_1tl + # CV admin 
    xpolity_1tl + anocracy_1tl + crgdppc_ln_1tl + cgdpg_1tl + cpop_ln_1tl + cpopg_1tl + incidencev412 + # CV country
    peaceyears_sc + peaceyears_sc_sq + peaceyears_sc_cb, # time dep
  data = cvdm_restricted_df)
coeftest(s3.m4,vcovHC(s3.m4,"HC0"))
# N obs
nobs(s3.m4)
# Dispersion parameter (theta) + associated SE
s3.m4$theta; s3.m4$SE.theta
# Loglik
logLik(s3.m4)
# aic
AIC(s3.m4)


# Model s4.m4: conditional fixed-effect logit
s4.m4 <- clogit(sc_inc_all ~ DIDP_1000pw + iterm1_pw + floodpw + iterm2_pw +# full overlap (t0) 
    DIDP_1000npw + iterm1_npw + floodnpw + iterm2_npw + # part. overlap (t0) 
    DIDP_1000near + iterm1_near + floodnear + iterm2_near + # near (t0) 
    DIDP_1000pw_1tl + iterm1_pw_1tl + floodpw_1tl + iterm2_pw_1tl +# full overlap (t-1)
    DIDP_1000npw_1tl + iterm1_npw_1tl + floodnpw_1tl + iterm2_npw_1tl + # part. overlap (t-1)
    DIDP_1000near_1tl + iterm1_near_1tl + floodnear_1tl + iterm2_near_1tl + # near (t-1)
    idp_ra + # conflict displacement variable
    flood_pexp + flood_pexp_sq + # count past events (all events)
    agdppc_ln + iapop_ln + sc_inc_all_1tl  + sc_cl_inc_1tl + # CV admin 
    xpolity_1tl + anocracy_1tl + crgdppc_ln_1tl + cgdpg_1tl + cpop_ln_1tl + cpopg_1tl + incidencev412 + # CV country
    peaceyears_sc + peaceyears_sc_sq + peaceyears_sc_cb + # time dep
    strata(admin1_code),
  data = cvdm_restricted_df)
summary(s4.m4)
# N obs
nrow(model.matrix(s4.m4))
# Loglik
logLik(s4.m4)
# aic
AIC(s4.m4)


## Model s5.m4: UCDP non-state conflict
s5.m4 <- bayesglm(ns_inc ~ DIDP_1000pw + iterm1_pw + floodpw + iterm2_pw +# full overlap (t0) 
    DIDP_1000npw + iterm1_npw + floodnpw + iterm2_npw + # part. overlap (t0) 
    DIDP_1000near + iterm1_near + floodnear + iterm2_near + # near (t0)
    DIDP_1000pw_1tl + iterm1_pw_1tl + floodpw_1tl + iterm2_pw_1tl +# full overlap (t-1)
    DIDP_1000npw_1tl + iterm1_npw_1tl + floodnpw_1tl + iterm2_npw_1tl + # part. overlap (t-1)
    DIDP_1000near_1tl + iterm1_near_1tl + floodnear_1tl + iterm2_near_1tl + # near (t-1)
    idp_ra + # conflict displacement variable
    flood_pexp + flood_pexp_sq + # count past events (all events)
    agdppc_ln + iapop_ln + ns_inc_1tl  + ns_cl_inc_1tl + # CV admin 
    xpolity_1tl + anocracy_1tl + crgdppc_ln_1tl + cgdpg_1tl + cpop_ln_1tl + cpopg_1tl + incidencev412 + # CV country
    peaceyears_ns + peaceyears_ns_sq + peaceyears_ns_cb, # time dep
  family = binomial(link="logit"),
  data = cvdm_restricted_df)
coeftest(s5.m4,vcovHC(s5.m4,"HC0"))
# N obs
nobs(s5.m4)
# Loglik
logLik(s5.m4)
# aic
AIC(s5.m4)

# Figure A.10: average predictive differences
apd_matrix_m4S5_0tl<-qi_inter(s5.m4,lag=F)
apd_matrix_m4S5_1tl<-qi_inter(s5.m4,lag=T)
par(ps=15,las=1,mar=c(0,0,0,2),oma=c(4, 12, 3, 0),mfcol=c(1,2))
# plot 1
plot.new();plot.window(xlim = limx,ylim=range(yat));abline(v=0,col="grey")
plotCI(apd_matrix_m4S5_1tl[,1],yat,
    li=apd_matrix_m4S5_1tl[,3],ui=apd_matrix_m4S5_1tl[,6],err='x',
    gap=0,xlab="",bty="n",
    ylab="",yaxt="n",pch=16,add=T);axis(1,cex.axis=.6)
axis(side=2, at=yat, labels=c(
    "Fully affected, Floods","Fully affected, Displ.",expression(italic("Difference")),
    "Partially affected, Floods","Partially affected, Displ.",expression(italic("Difference")),
    "Proximate, Floods","Proximate, Displ.",expression(italic("Difference"))
    ),cex.axis=.6,tick=FALSE)
mtext(expression(paste(Delta," P (Non-state conflict)")), 1, 2.5, cex=.6, outer=FALSE)
mtext("Lagged response (T-1)", 3, 1, cex=.6,at=0, outer=FALSE)
# plot 2
plot.new();plot.window(xlim = limx,ylim=range(yat));abline(v=0,col="grey")
plotCI(apd_matrix_m4S5_0tl[,1],yat,
    li=apd_matrix_m4S5_0tl[,3],ui=apd_matrix_m4S5_0tl[,6],err='x',
    gap=0,xlab="",bty="n",
    ylab="",yaxt="n",pch=16,add=T);axis(1,cex.axis=.6)
mtext(expression(paste(Delta," P (Non-state conflict)")), 1, 2.5, cex=.6, outer=FALSE)
mtext("Immediate response (T0)", 3, 1, cex=.6,at=0, outer=FALSE)


## Model s6.m4: UCDP civil war conflict
s6.m4 <- bayesglm(sb_inc ~ DIDP_1000pw + iterm1_pw + floodpw + iterm2_pw +# full overlap (t0)
    DIDP_1000npw + iterm1_npw + floodnpw + iterm2_npw + # part. overlap (t0) 
    DIDP_1000near + iterm1_near + floodnear + iterm2_near + # near (t0)
    DIDP_1000pw_1tl + iterm1_pw_1tl + floodpw_1tl + iterm2_pw_1tl +# full overlap (t-1)
    DIDP_1000npw_1tl + iterm1_npw_1tl + floodnpw_1tl + iterm2_npw_1tl + # part. overlap (t-1)
    DIDP_1000near_1tl + iterm1_near_1tl + floodnear_1tl + iterm2_near_1tl + # near (t-1)
    idp_ra + # conflict displacement variable
    flood_pexp + flood_pexp_sq + # count past events (all events)
    agdppc_ln + iapop_ln + sb_inc_1tl  + sb_cl_inc_1tl + # CV admin 
    xpolity_1tl + anocracy_1tl + crgdppc_ln_1tl + cgdpg_1tl + cpop_ln_1tl + cpopg_1tl + incidencev412 + # CV country
    peaceyears_sb + peaceyears_sb_sq + peaceyears_sb_cb, # time dep
  family = binomial(link="logit"),
  data = cvdm_restricted_df)
coeftest(s6.m4,vcovHC(s6.m4,"HC0"))
# N obs
nobs(s6.m4)
# Loglik
logLik(s6.m4)
# aic
AIC(s6.m4)

# Figure A.11: average predictive differences
apd_matrix_m4S6_0tl<-qi_inter(s6.m4,lag=F)
apd_matrix_m4S6_1tl<-qi_inter(s6.m4,lag=T)
par(ps=15,las=1,mar=c(0,0,0,2),oma=c(4, 12, 3, 0),mfcol=c(1,2))
# plot 1
plot.new();plot.window(xlim = limx,ylim=range(yat));abline(v=0,col="grey")
plotCI(apd_matrix_m4S6_1tl[,1],yat,
    li=apd_matrix_m4S6_1tl[,3],ui=apd_matrix_m4S6_1tl[,6],err='x',
    gap=0,xlab="",bty="n",
    ylab="",yaxt="n",pch=16,add=T);axis(1,cex.axis=.6)
axis(side=2, at=yat, labels=c(
    "Fully affected, Floods","Fully affected, Displ.",expression(italic("Difference")),
    "Partially affected, Floods","Partially affected, Displ.",expression(italic("Difference")),
    "Proximate, Floods","Proximate, Displ.",expression(italic("Difference"))
    ),cex.axis=.6,tick=FALSE)
mtext(expression(paste(Delta," P (Civil war)")), 1, 2.5, cex=.6, outer=FALSE)
mtext("Lagged response (T-1)", 3, 1, cex=.6,at=0, outer=FALSE)
# plot 2
plot.new();plot.window(xlim = limx,ylim=range(yat));abline(v=0,col="grey")
plotCI(apd_matrix_m4S6_0tl[,1],yat,
    li=apd_matrix_m4S6_0tl[,3],ui=apd_matrix_m4S6_0tl[,6],err='x',
    gap=0,xlab="",bty="n",
    ylab="",yaxt="n",pch=16,add=T);axis(1,cex.axis=.6)
mtext(expression(paste(Delta," P (Civil war)")), 1, 2.5, cex=.6, outer=FALSE)
mtext("Immediate response (T0)", 3, 1, cex=.6,at=0, outer=FALSE)